home *** CD-ROM | disk | FTP | other *** search
/ Space & Astronomy / Space and Astronomy (October 1993).iso / mac / programs / space / AA_51.ZIP / KEPLER.C < prev    next >
C/C++ Source or Header  |  1993-02-13  |  10KB  |  410 lines
  1.  
  2. /* Program to solve Keplerian orbit
  3.  * given orbital parameters and the time.
  4.  * Returns Heliocentric equatorial rectangular coordinates of
  5.  * the object.
  6.  *
  7.  * This program detects several cases of given orbital elements.
  8.  * If a program for perturbations is pointed to, it is called
  9.  * to calculate all the elements.
  10.  * If there is no program, then the mean longitude is calculated
  11.  * from the mean anomaly and daily motion.
  12.  * If the daily motion is not given, it is calculated
  13.  * by Kepler's law.
  14.  * If the eccentricity is given to be 1.0, it means that
  15.  * meandistance is really the perihelion distance, as in a comet
  16.  * specification, and the orbit is parabolic.
  17.  *
  18.  * Reference: Taff, L.G., "Celestial Mechanics, A Computational
  19.  * Guide for the Practitioner."  Wiley, 1985.
  20.  */
  21.  
  22. #include "kep.h"
  23.  
  24. #if !BandS
  25. extern struct orbit earth;    /* orbital elements of the earth */
  26. #endif
  27.  
  28. extern double eps, coseps, sineps; /* obliquity of ecliptic */
  29. double sinh(), cosh(), tanh();
  30. int embofs();
  31.  
  32. int kepler(J, e, rect, polar)
  33. double J, rect[], polar[];
  34. struct orbit *e;
  35. {
  36. double alat, E, M, W, temp;
  37. double epoch, inclination, ascnode, argperih;
  38. double meandistance, dailymotion, eccent, meananomaly;
  39. double r, coso, sino, cosa;
  40. int k;
  41.  
  42. /* Compute orbital elements if a program for doing so
  43.  * is supplied
  44.  */
  45. if( e->oelmnt )
  46.     {
  47.     k = (*(e->oelmnt) )(e,J);
  48.     if( k == -1 )
  49.         goto dobs;
  50.     }
  51. else if( e->celmnt )
  52.     { /* call B & S algorithm */
  53. dobs:
  54.     (*(e->celmnt) )( J, polar );
  55.     polar[0] = modtp( polar[0] );
  56.     E = polar[0]; /* longitude */
  57.     e->L = E;
  58.     W = polar[1]; /* latitude */
  59.     r = polar[2]; /* radius */
  60.     e->r = r;
  61.     e->epoch = J;
  62.     e->equinox = J;
  63.     goto kepdon;
  64.     }
  65.  
  66. /* Decant the parameters from the data structure
  67.  */
  68. epoch = e->epoch;
  69. inclination = e->i;
  70. ascnode = e->W * DTR;
  71. argperih = e->w;
  72. meandistance = e->a; /* semimajor axis */
  73. dailymotion = e->dm;
  74. eccent = e->ecc;
  75. meananomaly = e->M;
  76. /* Check for parabolic orbit. */
  77. if( eccent == 1.0 )
  78.     {
  79. /* meandistance = perihelion distance, q
  80.  * epoch = perihelion passage date
  81.  */
  82.     temp = meandistance * sqrt(meandistance);
  83.     W = (J - epoch ) * 0.0364911624 / temp;
  84. /* The constant above is 3 k / sqrt(2),
  85.  * k = Gaussian gravitational constant = 0.01720209895 . */
  86.     E = 0.0;
  87.     M = 1.0;
  88.     while( fabs(M) > 1.0e-11 )
  89.         {
  90.         temp = E * E;
  91.         temp = (2.0 * E * temp + W)/( 3.0 * (1.0 + temp));
  92.         M = temp - E;
  93.         if( temp != 0.0 )
  94.             M /= temp;
  95.         E = temp;
  96.         }
  97.     r = meandistance * (1.0 + E * E );
  98.     M = atan( E );
  99.     M = 2.0 * M;
  100.     alat = M + DTR*argperih;
  101.     goto parabcon;
  102.     }
  103. if( eccent > 1.0 )
  104.     {
  105. /* The equation of the hyperbola in polar coordinates r, theta
  106.  * is r = a(e^2 - 1)/(1 + e cos(theta))
  107.  * so the perihelion distance q = a(e-1),
  108.  * the "mean distance"  a = q/(e-1).
  109.  */
  110.     meandistance = meandistance/(eccent - 1.0);
  111.     temp = meandistance * sqrt(meandistance);
  112.     W = (J - epoch ) * 0.01720209895 / temp;
  113. /* solve M = -E + e sinh E */
  114.     E = W/(eccent - 1.0);
  115.     M = 1.0;
  116.     while( fabs(M) > 1.0e-11 )
  117.         {
  118.         M = -E + eccent * sinh(E) - W;
  119.         E += M/(1.0 - eccent * cosh(E));
  120.         }
  121.     r = meandistance * (-1.0 + eccent * cosh(E));
  122.     temp = (eccent + 1.0)/(eccent - 1.0);
  123.     M = sqrt(temp) * tanh( 0.5*E );
  124.     M = 2.0 * atan(M);    
  125.     alat = M + DTR*argperih;
  126.     goto parabcon;
  127.     }
  128. /* Calculate the daily motion, if it is not given.
  129.  */
  130. if( dailymotion == 0.0 )
  131.     {
  132. /* The constant is 180 k / pi, k = Gaussian gravitational constant.
  133.  * Assumes object in heliocentric orbit is massless.
  134.  */
  135.     dailymotion = 0.9856076686/(e->a*sqrt(e->a));
  136.     }
  137. dailymotion *= J - epoch;
  138. /* M is proportional to the area swept out by the radius
  139.  * vector of a circular orbit during the time between
  140.  * perihelion passage and Julian date J.
  141.  * It is the mean anomaly at time J.
  142.  */
  143. M = DTR*( meananomaly + dailymotion );
  144. M = modtp(M);
  145. /* If mean longitude was calculated, adjust it also
  146.  * for motion since epoch of elements.
  147.  */
  148. if( e->L )
  149.     {
  150.     e->L += dailymotion;
  151.     e->L = mod360( e->L );
  152.     }
  153.  
  154. /* By Kepler's second law, M must be equal to
  155.  * the area swept out in the same time by an
  156.  * elliptical orbit of same total area.
  157.  * Integrate the ellipse expressed in polar coordinates
  158.  *     r = a(1-e^2)/(1 + e cosW)
  159.  * with respect to the angle W to get an expression for the
  160.  * area swept out by the radius vector.  The area is given
  161.  * by the mean anomaly; the angle is solved numerically.
  162.  * 
  163.  * The answer is obtained in two steps.  We first solve
  164.  * Kepler's equation
  165.  *    M = E - eccent*sin(E)
  166.  * for the eccentric anomaly E.  Then there is a
  167.  * closed form solution for W in terms of E.
  168.  */
  169.  
  170. E = M; /* Initial guess is same as circular orbit. */
  171. temp = 1.0;
  172. do
  173.     {
  174. /* The approximate area swept out in the ellipse */
  175.     temp = E - eccent * sin(E)
  176. /* ...minus the area swept out in the circle */
  177.         - M;
  178. /* ...should be zero.  Use the derivative of the error
  179.  * to converge to solution by Newton's method.
  180.  */
  181.     E -= temp/(1.0 - eccent*cos(E));
  182.     }
  183. while( fabs(temp) > 1.0e-11 );
  184.  
  185. /* The exact formula for the area in the ellipse is
  186.  *    2.0*atan(c2*tan(0.5*W)) - c1*eccent*sin(W)/(1+e*cos(W))
  187.  * where
  188.  *    c1 = sqrt( 1.0 - eccent*eccent )
  189.  *    c2 = sqrt( (1.0-eccent)/(1.0+eccent) ).
  190.  * Substituting the following value of W
  191.  * yields the exact solution.
  192.  */
  193. temp = sqrt( (1.0+eccent)/(1.0-eccent) );
  194. W = 2.0 * atan( temp * tan(0.5*E) );
  195.  
  196. /* The true anomaly.
  197.  */
  198. W = modtp(W);
  199.  
  200. meananomaly *= DTR;
  201. /* Orbital longitude measured from node
  202.  * (argument of latitude)
  203.  */
  204. if( e->L )
  205.     alat = (e->L)*DTR + W - meananomaly - ascnode;
  206. else
  207.     alat = W + DTR*argperih; /* mean longitude not given */
  208.  
  209. /* From the equation of the ellipse, get the
  210.  * radius from central focus to the object.
  211.  */
  212. r = meandistance*(1.0-eccent*eccent)/(1.0+eccent*cos(W));
  213.  
  214. parabcon:
  215.  
  216. /* The heliocentric ecliptic longitude of the object
  217.  * is given by
  218.  *   tan( longitude - ascnode )  =  cos( inclination ) * tan( alat ).
  219.  */
  220. coso = cos( alat );
  221. sino = sin( alat );
  222. inclination *= DTR;
  223. W = sino * cos( inclination );
  224. E = zatan2( coso, W ) + ascnode;
  225.  
  226. /* The ecliptic latitude of the object
  227.  */
  228. W = sino * sin( inclination );
  229. W = asin(W);
  230.  
  231. /* Apply perturbations, if a program is supplied.
  232.  */
  233. if( e->celmnt )
  234.     {
  235.     e->L = E;
  236.     e->r = r;
  237.     (*(e->celmnt) )(e);
  238.     E = e->L;
  239.     r = e->r;
  240.     W += e->plat;
  241.     }
  242.  
  243. /* If earth, Adjust from earth-moon barycenter to earth
  244.  * by AA page E2
  245.  * unless orbital elements are calculated by formula.
  246.  * (The Meeus perturbation formulas include this term for the moon.)
  247.  */
  248. #if !BandS
  249. if( (e == &earth) && (e->oelmnt == 0) )
  250.     {
  251.     embofs( J, &E, &W, &r ); /* see below */
  252.     }
  253. #endif
  254.  
  255. /* Output the polar cooordinates
  256.  */
  257. polar[0] = E; /* longitude */
  258. polar[1] = W; /* latitude */
  259. polar[2] = r; /* radius */
  260.  
  261.  
  262. kepdon:
  263.  
  264.  
  265. /* Convert to rectangular coordinates,
  266.  * using the perturbed latitude.
  267.  */
  268. rect[2] = r * sin(W);
  269. cosa = cos(W);
  270. rect[1] = r * cosa * sin(E);
  271. rect[0] = r * cosa * cos(E);
  272.  
  273. /* Convert from heliocentric ecliptic rectangular
  274.  * to heliocentric equatorial rectangular coordinates
  275.  * by rotating eps radians about the x axis.
  276.  */
  277. epsiln( e->equinox );
  278. W = coseps*rect[1] - sineps*rect[2];
  279. M = sineps*rect[1] + coseps*rect[2];
  280. rect[1] = W;
  281. rect[2] = M;
  282.  
  283. /* Precess the position
  284.  * to ecliptic and equinox of J2000.0
  285.  * if not already there.
  286.  */
  287. precess( rect, e->equinox, 1 );
  288. return(0);
  289. }
  290.  
  291.  
  292. #if !BandS
  293. /* Adjust position from Earth-Moon barycenter to Earth
  294.  *
  295.  * J = Julian day number
  296.  * E = Earth's ecliptic longitude (radians)
  297.  * W = Earth's ecliptic latitude (radians)
  298.  * r = Earth's distance to the Sun (au)
  299.  */
  300. int embofs( J, E, W, r )
  301. double J;
  302. double *E, *W, *r;
  303. {
  304. double T, M, a, L, B, p;
  305. double smp, cmp, s2mp, c2mp, s2d, c2d, sf, cf;
  306. double s2f, sx, cx, xyz[3];
  307.  
  308. /* Short series for position of the Moon
  309.  */
  310. T = (J-J1900)/36525.0;
  311. /* Mean anomaly of moon (MP) */
  312. a = ((1.44e-5*T + 0.009192)*T + 477198.8491)*T + 296.104608;
  313. a *= DTR;
  314. smp = sin(a);
  315. cmp = cos(a);
  316. s2mp = 2.0*smp*cmp;        /* sin(2MP) */
  317. c2mp = cmp*cmp - smp*smp;    /* cos(2MP) */
  318. /* Mean elongation of moon (D) */
  319. a = ((1.9e-6*T - 0.001436)*T + 445267.1142)*T + 350.737486;
  320. a  = 2.0 * DTR * a;
  321. s2d = sin(a);
  322. c2d = cos(a);
  323. /* Mean distance of moon from its ascending node (F) */
  324. a = (( -3.e-7*T - 0.003211)*T + 483202.0251)*T + 11.250889;
  325. a  *= DTR;
  326. sf = sin(a);
  327. cf = cos(a);
  328. s2f = 2.0*sf*cf;    /* sin(2F) */
  329. sx = s2d*cmp - c2d*smp;    /* sin(2D - MP) */
  330. cx = c2d*cmp + s2d*smp;    /* cos(2D - MP) */
  331. /* Mean longitude of moon (LP) */
  332. L = ((1.9e-6*T - 0.001133)*T + 481267.8831)*T + 270.434164;
  333. /* Mean anomaly of sun (M) */
  334. M = (( -3.3e-6*T - 1.50e-4)*T + 35999.0498)*T + 358.475833;
  335. /* Ecliptic longitude of the moon */
  336. L =    L
  337.     + 6.288750*smp
  338.     + 1.274018*sx
  339.     + 0.658309*s2d
  340.     + 0.213616*s2mp
  341.     - 0.185596*sin( DTR * M )
  342.     - 0.114336*s2f;
  343. /* Ecliptic latitude of the moon */
  344. a = smp*cf;
  345. sx = cmp*sf;
  346. B =      5.128189*sf
  347.     + 0.280606*(a+sx)        /* sin(MP+F) */
  348.     + 0.277693*(a-sx)        /* sin(MP-F) */
  349.     + 0.173238*(s2d*cf - c2d*sf);    /* sin(2D-F) */
  350. B *= DTR;
  351. /* Parallax of the moon */
  352. p =     0.950724
  353.     +0.051818*cmp
  354.     +0.009531*cx
  355.     +0.007843*c2d
  356.     +0.002824*c2mp;
  357. p *= DTR;
  358. /* Elongation of Moon from Sun
  359.  */
  360. L = DTR * mod360(L) - (*E - PI);
  361.  
  362. #if !OFDATE
  363. /* The Moon's position has to be precessed to J2000 */
  364.  
  365. /* Distance in au */
  366. a = 4.263523e-5/sin(p);
  367.  
  368. /* Convert to rectangular ecliptic coordinates */
  369. xyz[2] = a * sin(B);
  370. cx = cos(B);
  371. xyz[1] = a * cx * sin(L);
  372. xyz[0] = a * cx * cos(L);
  373.  
  374. /* Convert to equatorial */
  375. epsiln( J );
  376. sf = coseps*xyz[1] - sineps*xyz[2];
  377. cf = sineps*xyz[1] + coseps*xyz[2];
  378. xyz[1] = sf;
  379. xyz[2] = cf;
  380.  
  381. /* Precess the position
  382.  * to ecliptic and equinox of J2000.0
  383.  */
  384. precess( xyz, J, 1 );
  385.  
  386. /* Rotate back to ecliptic coordinates */
  387. sf = coseps*xyz[1] + sineps*xyz[2];
  388. cf = -sineps*xyz[1] + coseps*xyz[2];
  389. xyz[1] = sf;
  390. xyz[2] = cf;
  391.  
  392. /* Convert to polar coordinates */
  393. L = zatan2( xyz[0], xyz[1] );
  394. B = asin( xyz[2]/a );
  395. #endif
  396.  
  397. /* Distance from Earth to Earth-Moon barycenter
  398. */
  399. a = 5.18041e-7 / p;
  400.  
  401. /* Adjust the coordinates of the Earth
  402.  */
  403. sx = a / *r;        /* chord = R * dTheta */
  404. *E += sx * sin(L);
  405. *W -= sx * B;
  406. *r += a * cos(L);
  407. return(0);
  408. }
  409. #endif
  410.